Quantum statistics and the momentum distribution of liquid para-hydrogen 
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Extensive Monte Carlo simulations of bulk liquid para-hydrogen at a temperature T=16.5 K have 
been carried out using the continuous-space Worm Algorithm. Results for the momentum distribu- 
tion, as well as for the kinetic energy per particle and the pair correlation function are reported. The 
static equilibrium thermodynamic properties of this system can be generally computed by assum- 
ing that molecules are distinguishable. However, the one-body density matrix (and the associated 
momentum distribution) are affected by particle indistinguishability and quantum statistics, to an 
extent that lends itself to experimental observation. Comparisons with available experimental data 
and other theoretical and numerical calculations are offered. 

PACS numbers: 67.10.Hk, 61.20.Ja, 67.63.Cd, 67. 
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I. INTRODUCTION 

The investigation of the condensed phases of molecular 
hydrogen is aimed at understanding and characterizing 
quantitatively quantum effects, both in the liquid and 
in the solid phases. In many respects, liquid hydrogen 
displays a physical behavior that interpolates between 
that of a classical liquid, and that of supcrfiuid helium, 
i.e., the most quantum-mechanical of the simple fluids. 
The most important manifestation of quantum mechan- 
ics in a fluid of hydrogen molecules, is the the fact that 
the kinetic energy per particle 1 significantly exceeds the 
classically predicted^ value 3T/2. This is essentially a 
consequence of zero-point motion, as the excess kinetic 
energy arises from the confinement that each molecule ex- 
periences inside the instantaneous "cage" of surrounding 
particles, owing to the hard core of the intcr-molccular 
interaction. 

Hardly any signature of quantum statistics (Bose in 
this both para- and ort/io-hydrogen molecules 

have integer spin) can be detected on the equilibrium 
thermodynamic properties of liquid hydrogen, as parti- 
cle exchanges are exceedingly rare. For a simple quanti- 
tative estimate, one may consider the lightest hydrogen 
isotope, i.e., para- hydrogen (P-H2); under the pressure 
of its own vapor, liquid P-H2 crystallizes at a tempera- 
ture T=13.8 K. At temperatures slightly above that, in 
the liquid phase, the thermal wavelength of a P-H2 
molecule is of the order of 1 A, significantly smaller than 
the mean intermolecular distance d ~ 3.5 A. Thus, prob- 
ability "clouds" associated to different molecules do not 
overlap significantly, and exchanges are suppressed 3 - at 
least as ~ exp[— d 2 /2A^], which is less than 10~ 3 (for 
comparison, it is close to 0.3 in 4 He at T=l K). 

The presence of the above-mentioned hard core repul- 
sion at distances less than ~ 2.2 A in the interaction 
between two molecules, suppresses quantum exchanges 
even further^ Consequently, liquid hydrogen has been 
modeled as a system of distinguishable quantum parti- 
cles ("Boltzmannons"), e.g., in all numerical simulation 
work to date Indeed, it is precisely the impos- 

sibility of stabilizing a phase of this system where quan- 
tum exchanges may be important jl! that has so far ham- 
pered the experimental observation of superfluidity in 



bulk para-hydrogen. 

However, effects of particle indistinguishability and 
quantum statistics ought to be detectable in the momen- 
tum distribution, defined as 



n(k) = (a^ak) 



(1) 



where (...) stands for thermal expectation value, and d^ 
(dk) is the creation (annihilation) operator for a particle 
of momentum k. The momentum distribution is related, 
via a Fourier transformation, to the one-body density ma- 
trix 



(2) 



where ip and are Bose field operators. Within the path 
integral formulation of quantum statistical mechanics, 3 
n(r,r') describes the distribution of relative positions of 
the two dangling ends of a single-particle path that has 
been "cut open". 

Due to quantum exchanges (even as rare as they are 
in liquid hydrogen), and the ensuing entanglement of 
single-particle paths, n(|r — r'|) extends out to signifi- 
cantly greater distances than it would if particles were 
truly distinguishable (in which case n(r) is very nearly 
Gaussian) . This in turn results in a transfer of weight of 
n(k) toward lower momenta, which can be interpreted as 
a sign of the incipient Bose-Einstein Condensation that 
liquid para-hydrogen would undergo, at a temperature 
around 5 K, were it to escape crystallization. It has been 
suggested^ 2 - that this effect ought to be experimentally 
measurable. 

The bulk of the available experimental information on 
the structure and dynamics of liquid hydrogen comes 
from Deep Inelastic Neutron Scattering (DINS) measure- 
ments, which offer (in some cases direct) access to quan- 
tities such as the single-particle kinetic energy, momen- 
tum distribution and excitation spectrum^ Such exper- 
iments have been pursued by various groups, over the 
past two decadesJ^±£i±&±!i±&±&2ii22 Qn the theoreti _ 

cal side, besides the above-mentioned numerical simula- 
tions, typically based on Quantum Monte Carlo (QMC) 
techniques, an analytical approach known as Correlated 
Density Matrix (CDM) theory^^ has been shown to 
provide not only qualitative insight, but also quantita- 
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tively reliable predictions, at least for structural prop- 
erties (an additional check on the predictions based on 
CDM is furnished here). 

Agreement between theoretical predictions and exper- 
imental data has not been of the quality that one would 
expect, given the availability of fairly well-established 
simulation methods and reasonably quantitative micro- 
scopic models. For example, the difference between the 
kinetic energy per particle in liquid p-fi2 at T=16.5 K, 
in a range of pressure between 1 and 80 bars, computed 
theoretically 4, and inferred from the most recent experi- 
mental measurements 1 ^ is approximately 10% (a few K), 
which seems large, considering that the same compari- 
son yields much better agreement for the more quantal 
helium liquid. — It is unclear whether such a discrepancy 
originates within the microscopic model, or the computa- 
tional methodologies, or may lie instead with the analysis 
of the experimental data. No results of any first princi- 
ples microscopic calculation of the momentum distribu- 
tion in liquid p-H 2 have been reported so far. 

The purpose of this paper is, on the one hand, to 
provide independent theoretical estimates for the kinetic 
energy per particle, as well as for the pair correlation 
function at the same conditions of temperature and pres- 
sure considered in Refs. |4j and [12, thereby allowing for 
an extended comparison between different calculations. 
More importantly, an unbiased theoretical estimate of 
the single-particle momentum distribution is furnished 
here, enabling a direct and cogent comparison between 
theory and experiment. The computational tool uti- 
lized here is numerical; specifically, use is made of the 
continuous-space Worm Algorithm (WA), which allows 
for a direct calculation of the one-particle density matrix, 
connected to the momentum distribution via a straight- 
forward Fourier transformation. 

Our numerical estimates for the kinetic energy are for 
the most part in quantitative agreement with those of 
other works, and thus retain the same level of disagree- 
ment with experimental data already reported by others. 
Quantum effects are clearly seen in the computed one- 
body density matrix, which can be directly associated 
to quantum exchanges, explicitly allowed in our calcu- 
lation (i.e., no assumption of distinguishability is made 
here). In the next section, the microscopic model under- 
lying the calculation as well as the basic features of the 
methodology utilized are reviewed. In Se. IIIII a thor- 
ough illustration of the results obtained in this work is 
provided; finally, a general discussion is offered, and con- 
clusions outlined, in Sec. IIVI 



II. MODEL 

Consistently with all previous theoretical studies, our 
system of interest is modeled as an ensemble of N p- 
H2 molecules, regarded as point particles of spin zero, 
enclosed in a cubic vessel of volume fl with periodic 
boundary conditions. The quantum-mechanical many- 



body Hamiltonian is the following: 

N 

i=l i<j 

Here, A=12.031 KA 2 for p-U 2 , while V is the potential 
describing the interaction between two molecules, only 
depending on their relative distance. Most of the results 
presented here were obtained using the Silvera-Goldman 
model potential. 24 This is not the only potential that has 
been used in previous work, but it is arguably the most 
commonly adopted. It also been showr* 2 ^ to afford an 
accurate quantitative description of the thermodynam- 
ics of the solid phase of p-H2- For the sake of compari- 
son, however, we have also performed a calculation using 
the Buck 2 - i 27 ' 28 pair potential. Naturally, in principle a 
more accurate model would go beyond the simple pair 
decomposition, including, for instance, also interactions 
among triplets; however, published numerical work (e.g., 
on helium) has given strong indications that three-body 
corrections, while significantly affecting the estimation of 
the pressure, have a relatively small effect on the struc- 
ture and dynamics of the system, of interest here . 29 ' 30 

The thermodynamic properties of the system, as mod- 
eled by the many-body Hamiltonian ([3]) , were studied by 
means of numerical simulations, based on the continuous- 
space Worm Algorithm (WA) ,22-22 This (Monte Carlo) 
methodology, recently introduced, has several advantages 
over Path Integral Monte Carlo (PIMC), the technique 
utilized in most previous many-body numerical calcula- 
tions for liquid P-H2. The most important, in this con- 
text, is the fact that it is formulated in an extended con- 
figuration space, including an open world line. This al- 
lows one to evaluate, besides all usual thermodynamic 
observables, also off-diagonal single-particle correlation 
functions not accessible in conventional PIMC, at no 
added computational cost. One of these correlation func- 
tions is the single-particle Matsubara Green function 

g(r,r',T) = (f[^(v,r) ^(r',0)]} (4) 

Here, ip and ijp are time-dependent field operators, —/3 < 
t < f3, with (3 = l/T is commonly referred to as "imag- 
inary time" , and T is the time ordering operator.— In 
the limit r — > 0, the Matsubara Green function reduces 
to the one-particle density matrix n(r, r'), defined in Eq. 
©■ 

For a translationally invariant system, it is n(r',r) = 
ri(r' — r), and in three dimensions the momentum distri- 
bution n(k) is related to n(r) through 

n(r) = J d 3 k n(k) e* kr (5) 

For a system that is also isotropic, like a liquid, it is 
n(r) = ri(r). Hence, on inverting the above relationship 
one obtains 

4tt f°° 

nCk) = n(k) = —— I dr r sinfcr n(r) (6) 
k Jo 
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For convenience, the following normalization is imposed 
on n(k): 

^3 J d?k n(k) = 1 (7) 

which fixes to unity the value of the one-body density 
matrix at r = 0. The average kinetic energy per particle 
Ek is connected to the momentum distribution through 

1 r T, 2 h 2 

Ek = TTv d*k^n(k) (8) 
(27rj J J 2m 

Using ©, and ©, it is straightforward to show that, 
in the limit r — ■> 0, it is 

n(r) « 1 - §r 2 (9) 

which provides a useful consistency check on the com- 
puted Ek and n(r). We obtain n(k) through (J6]), using 
the numerically computed n(r). 

The reader is referred to Ref . [32] for a thorough descrip- 
tion of the continuous-space WA. The specific implemen- 
tation utilized in this project is canonical, i.e., we keep 
the number N of particles fixed, 34 Other technical as- 
pects of the calculations are common to any other QMC 
simulation scheme. Results were obtained by simulating 
systems comprising two different numbers of particles, 
namely 7V=96 and iV=256. Within the statistical errors 
of the calculations, no significant size dependence can be 
detected in any of the quantities considered here, except 
for the kinetic energy, for which a systematic difference 
of approximately 0.1 K exists for the two system sizes 
(the larger system yielding the greater value of Ek). 
The usual fourth-order high-temperature propagator uti- 
lized in all previous studies based on the WA was adopted 
here; convergence of the estimates was observed for a 
time step r « 1.89 x 1CP 3 K _1 , which corresponds to 
P=32 imaginary time "slices" at the single temperature 
considered in this work. Because the computational cost 
was negligible, all estimates reported here were obtained 
using twice as small a time step, in order to be on the 
safe side. 



III. RESULTS 

The thermodynamic conditions considered in this work 
are the same as in Ref. 0. Specifically, the tempera- 
ture is fixed at T=16.5 K, and three different densities 
are considered, namely p=0.02235 A" 3 , 0.02372 A" 3 and 
0.02413 A~ 3 , corresponding to pressures of 1, 40 and 80 
bars respectively. 

A. Structure 

Static correlation functions are yielded directly by the 
simulation. The results obtained here are largely con- 
sistent with those of previous MC calculations**^, and 



shall therefore not be discussed any further. The results 
at the lowest and highest densities considered here are 
shown in Fig. [T] The basic physical features of these 
functions have already been extensively described in the 
literature. Differences between different calculations, as 
well as between theory and experiment, if any, are not 
easily spotted^ 




2 4 6 8 10 
r(A) 

FIG. 1: (Color online). Pair correlation function g(r) in liquid 
p-H 2 at T=16.5 K and p=0.02235 A" 3 (1 bar, solid line) and 
p—0. 02413 A -3 (80 bars, dashed line). Statistical errors are 
not visible, on the scale used here. 



B. Kinetic Energy 



Density (A~ 3 ) 0.02235 0.02373 0.02413 

Silvera-Goldman 61.78(2) 64.99(5) 67.93(4) 

Buck 62.61(3) 

PIMC (Ref. 4) 62.7(4) 65.5(5) 69.6(5) 

PIMC (Ref. 19) 61.4(1) 

Expt. (Ref. 18) 60.3(6) 

Expt. (Ref. .12) 67.8(3) 73.5(4) 77.5(4) 



TABLE I: Kinetic energy per p-IU molecule Ek (in K) at 
T=16.5 K, computed in this work (first two rows) and in 
Refs. lil and fTsL as well as recently determined experimentally 
(Refs. [IH an d 03) ■ Statistical errors (in parentheses) are 
on the last digit. Theoretical estimates are obtained using 
the Silvera-Goldmann potential, except for that in the second 
row, which is based on the Buck potential. Experimental 
value from Ref. [l|| is obtained as an interpolation of data at 
T=15.7 K (p=0. 02252 A" 3 ) and T=17.2 K (p=0. 02210 A" 3 ). 

Table U summarizes the kinetic energy per P-H2 
molecule computed in this work; results quoted are for a 
system of iV=256 molecules. For comparison, results of 
PIMC calculations published in Refs. |j and KL9l are also 
shown, together with experimental data from Refs. [H, 
as well as the more recent ones from Ref. It should 
be noted that the kinetic energy in QMC is not computed 
through the momentum distribution, but rather using a 
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direct estimator. 

The comparison of theoretical estimates is not completely 
satisfactory, even after statistical and possible systematic 
uncertainties are properly taken into account. At the 
lowest density, the result offered in Ref. HH is very close 
to the one obtained here, and it is almost certain that 
the small difference is attributable to a relatively small 
number of time slices (P=64) used in Ref. where the 
primitive approximation for the high-temperature prop- 
agator was utilized. While such a number is sufficient to 
obtain reasonably accurate estimate of structural prop- 
erties (e.g., the pair correlation function), it was found in 
this work that at least twice as many slices are needed, 
in order to achieve convergence of kinetic energy esti- 
mates to the precision quoted here, if the primitive ap- 
proximation is used (it is worth repeating that all results 
presented here are obtained using a fourth-order propa- 
gator). 

On the other hand, the discrepancy between our esti- 
mates and those of Ref. H is definitely outside statistical 
errors at the highest density, and possibly at the lowest 
one as well, in spite of the relatively large errors quoted 
therein; in general, the estimates offered in Ref. 3 ap- 
pear to be above ours, by almost 1.5-2 K at the highest 
density. Details of the calculation carried out in Ref. U 
are not available to us at this time, and therefore it is 
unclear what the source of such a disagreement might 
be, given that the model potential utilized is the same, 
and differences in system sizes seem unlikely to generate 
a disaccord of this magnitude. 

In order to carry out a comparison of the numerical re- 
sults obtained here with those obtained by other meth- 
ods, we have also performed a simulation at a density 
/9=0.021 A -3 for a temperature T=16 K. For this ther- 
modynamic condition, an estimate of the kinetic energy 
per particle of 58.6 K was provided in Ref. [H, based on 
CDM theory Our result is 57.4 ± 0.1 K, which seems in 
reasonable agreement. 

A much more significant discrepancy exists between 
the theoretical and some of the available experimental 
estimates for the kinetic energy, namely those reported 
in Ref. 0, all of which are above the theoretical ones, by 
approximately 6 K at the lowest density, up to almost 9 K 
at the highest. Such a disagreement is puzzling, consid- 
ering the quantitative agreement between theoretically 
computed and experimentally measured kinetic energy 
per particle in liquid helium, 23 in which deviations from 
the classical behavior are much more pronounced than 
in liquid hydrogen. In light of the substantial closeness 
(if not downright agreement) of the various independent 
theoretical estimates, one may look at the intermolecular 
potential as one of the possible sources of the discrep- 
ancy with experiment. The Silvera-Goldman potential 
has been shown to provide a rather accurate description 
of the crystalline phase of p-H.2 at low temperature^ and 
therefore it seems unlikely that it would not be at least as 
adequate in the liquid phase, in which the consequences 
of the pair-wise spherical approximation should be even 
less important than in the crystal. 



carried out in this work based on a different model po- 
tential as well, namely the Buck potential, at the lowest 
density considered here; the value is shown in Table [U 
The Buck potential gives a slightly higher kinetic energy, 
by about 1 K, a considerably smaller difference than that 
between theory and experiment. It should be noted that 
both these potentials yield a theoretical value of the ki- 
netic energy per particle in the solid phase within ~ 1-2 
K of the experimentally determined onei 18 ' 20 It seems 
therefore unlikely that the choice of pair potential may 
account for the observed difference between the theoret- 
ical estimates and the experimental data of Ref. [l^- 

Another possibility is that the finite size of the sim- 
ulated system may result in an underestimation of the 
kinetic energy. However, as mentioned above calcula- 
tions performed in this work on a system of A=96 par- 
ticles, yield essentially the same estimate obtained with 
7V=256, the difference being of the order of 0.1 K at the 
most. Therefore, it seems safe to exclude the possibility 
of a significant size dependence. Based on these consider- 
ations, we surmise that the discrepancy likely originates 
with some of the assumptions underlying the analysis of 
the experimental data carried out in Ref. [J, which may 
have to be reconsidered. On this point, it is worth not- 
ing the substantial disaccord between the experimental 
estimates for the same quantity provided in Refs. [TU and 
Il8l (about eight times the combined uncertainties quoted 
by the two groups). The value given in Ref. [H while 
not in perfect agreement with theoretical estimates, is 
nonetheless much closer to them than that furnished in 
Ref. [U. 



Momentum Distribution 




FIG. 2: (Color online) Theoretically computed one-body den- 
sity matrix n(r) in liquid p-Hb at T=16.5 K and p—0. 02235 
A~ 3 (circles). When not shown, statistical errors are smaller 
than symbol sizes. Straight line is an exponential fit to the 
curve for r > 6 A. Boxes show the result for n(r) for distin- 
guishable quantum particles; it is essentially identical with the 
one for Bosons for r less than ~ 4 A, but deviates significantly 
from it for larger distances, as a result of quantum exchanges. 



For comparison purposes, the same calculation was 



Fig. [5] shows the computed (spherically averaged) one- 



5 



body density matrix n(r) for liquid P-H2 at the lowest 
density considered here, namely p—0. 02235 A~ 3 . Devi- 
ations from the classical (Gaussian) behavior are clear 
(at the same time, the logarithmic scale gives an idea 
of their magnitude). In particular, n(r) extends beyond 
the first coordination shell, whose radius is ~ 4 A, i.e., a 
molecule has a nonzero probability amplitude of exiting 
the "cage" formed by the surrounding molecules. This 
is due to quantum-mechanical exchanges, which have in- 
deed been observed in all of the simulation carried out 
in this work (in the so-called G-sector - see Ref. [U for 
details), and which allow the two dangling ends of the 
open path to drift further away from each other than if 
particles were distinguishable, by virtue of entanglement 
among different paths. 

In order to establish this fact more quantitatively, we 
have performed a simulation in which particles were as- 
sumed truly distinguishable. In technical terms, this 
means that the so-called "swap" move^ which allows 
for entanglement of single-particle paths in the presence 
of a single open world line, is inhibited. The result that 
one obtains in this case for n(r), also shown in Fig. [21 is 
virtually identical with the one yielded by the calculation 
in which exchanges are included, for distances less than 
~ 4 A; for greater distances, on the other hand, it differs 
significantly, as the one-body density matrix continues to 
decay monotonically. On the scale of Fig. [2J no signal 
appears for r > 5 A, in such a simulation. This estab- 
lishes that the structure of n(r) above 4 A is a genuine 
consequence of quantum (Bose) statistics. 

From the numerically computed n(r), one can obtain 
the experimentally observable momentum distribution 
n(k) based on ©, by means of a straightforward nu- 
merical integration. Obtaining good statistics for values 
of n(r) at distances greater than ~ 6 A becomes quickly 
impractical, due to the rapid decay of the function, which 
requires an exceedingly long simulation time in order to 
achieve a meaningfully small statistical uncertainty. This 
could in principle be an issue, when trying to evaluate ([6]). 
However, the contribution to n(k) coming from distances 
greater than 6 A can be easily estimated by fitting the 
portion of the curve for r > 6 A, for example to an expo- 
nentially decaying function (as shown in Fig. [J). As it 
turns out, such contribution (which is greatest at k—0) 
is barely worth 0.1% of the total value of n(k), regardless 
on the particular fitting function that one chooses, and 
therefore an accurate evaluation of © can be obtained 
by integrating only up to r=6 A. 

Fig. [3] shows the momentum distribution n(k) result- 
ing from the numerical integration of ([6]) based on the 
data for n(r) shown in Fig. [5] Although the system stud- 
ied is not a classical liquid, the overall shape of n(k) re- 
mains close to a Gaussian^ albeit one corresponding to 
a different, effective temperature, given by 2Ek/3 (shown 
for comparison in Fig. [5]). Such a Gaussian is essentially 
indistinguishable from the momentum distribution that 
one obtains on Fourier transforming the one-body den- 
sity matrix computed for distinguishable quantum par- 
ticles. There is some strength transferred to both lower 
and higher momenta, compared to what one would find 
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FIG. 3: (Color online). Momentum distribution n(k) in liq- 
uid p-H 2 at T=16.5 K and p=0.02235 A -3 . Solid line is the 
estimate obtained by numerical integration of Eq. [6] using 
the results for n(r) shown in Fig. [2] as explained in the text. 
Statistical and systematic errors are too small to show on the 
scale of the figure. Dashed line is a Gaussian function whose 
width is chosen so as to yield the kinetic energy per particle 
computed in the simulation, namely 61.78±0.02 K. 
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FIG. 4: (Color online). Momentum distribution n(k) in liquid 
p-H 2 at T=16.5 K and p=0. 02235 A -3 (1 bar, solid line), 
p=0.02327 A" 3 (40 bars, dashed line), and p=0. 02413 A~ 3 
(80 bars, dotted line). Statistical errors are not visible, on 
the scale used here. 



if the momentum distribution were indeed such a " renor- 
malized" Gaussian; more quantitatively, numerical inte- 
gration of Eq. ijHJ) up to momentum fc =4 A -1 , using 
the data for n(k) shown in Fig. [3[ yields approximately 
85% of the overall kinetic energy (as opposed to 92.8% 
if n(k) were a Gaussian), the rest coming from momenta 
higher than k a . This underscores the delicacy of extract- 
ing the single-particle kinetic energy from the experimen- 
tally measured momentum distribution, as an accurate 
determination of the tail is required, a fact which might 
help account some of the differences reported in the lit- 
erature between the numerically computed and experi- 
mentally determined kinetic energy^ 

The difference between the computed momentum dis- 
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FIG. 5: (Color online). Kinetic energy distribution k 2 n(k) 
in liquid p-H 2 at T=16.5 K and p=0. 02235 A -3 (1 bar, solid 
line). Statistical errors are not visible, on the scale used here. 
The same units are used as in Ref. \lj 

tribution and the model Gaussian, which can be at- 
tributed entirely to quantum statistics, is altogether 
rather small, and is most noticeable for momenta be- 
tween 1.5 and 2.5 A -1 , corresponding to interparticle dis- 
tances between 2.5 and 4 A. This conclusion was already 
stated in Ref. H3 , based on an analysis of experimental 
data for n(k). In Fig. [5] the kinetic energy distribution 
is shown, namely the quantity k 2 n(k), for which exper- 
imental data are reported in Ref. [l2| (the same units 
utilized therein are used). Within the statistical uncer- 
tainties of the calculation carried out here, there appears 
to be broad agreement between the results obtained in 
this work and the experimental data, but with some no- 
ticeable differences. In particular, the computed k 2 n(k) 
attains its maximum in correspondence of the momen- 
tum k = 1.735±0.005 A~\ in excellent agreement with 
the value 1.736 A~ x reported in Ref. [l^- On the other 
hand, around k w 4 A -1 the curve calculated here falls 
below the experimental one, approximately by a factor 
2.5. This is consistent with the fact that the experimen- 
tally determined kinetic energy is above that obtained 
here by approximately 10%. 

The momentum distribution at higher density (pres- 
sure) does not qualitatively change, with respect to that 
shown in Fig. [2] This is shown in Fig. [1J where the 
momentum distribution computed at the three different 
densities considered here is displayed. The most impor- 
tant change that occurs on raising the pressure is the loss 
of weight at low momenta, as the system becomes in- 
creasingly classical. Correspondingly, the tail of the one- 
body density matrix at long distances is suppressed, as 
exchanges are rarer and rarer, even if one single-particle 



world line is open. 



IV. CONCLUSIONS 

We have carried out first principle calculations of 
the momentum distribution of liquid para-hydrogen at 
T=16.5 K, at three different densities corresponding to 
pressures ranging between 1 and 80 bars. The results 
for the one-particle density matrix show that quantum- 
mechanical exchanges result in a longer tail of the one- 
body density matrix than one would predict based on 
distinguishability of molecules. Consequently, this has 
an effect on the momentum distribution, which features 
clear, measurable deviations from a Gaussian. These con- 
clusions are in broad qualitative and quantitative agree- 
ment with considerations made in Ref. [12} On increas- 
ing the pressure, exchanges are suppressed, and the most 
significant change in the momentum distribution occurs 
precisely at low k. 

For the thermodynamic conditions explored here, ex- 
perimental measurements have been recently carried out. 
The agreement between theoretical estimates and exper- 
imental data continue to be less than satisfactory, cer- 
tainly of much lesser quality than that found for liquid 
and solid helium. The deviation between theoretical and 
experimental data appears to be systematic. At the time 
of this writing, it is unclear to us where the origin of the 
disagreement lies, but perhaps an independent check of 
the analysis of the data of Ref. [H] is in order. 

As mentioned above, the discrepancy could also be at- 
tributed to an inaccurate determination of the tail of the 
momentum distribution. However, another aspect worth 
revisiting is the effect of multiple inelastic scattering that 
a neutron suffers in a liquid. An example of the impact 
of multiple scattering is discussed in Ref. H3, where it is 
shown that it can indeed lead to an overestimation of the 
value of the center-of-mass kinetic energy per particle. 
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